Test of the semischematic model for a liquid of linear molecules 
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We apply to a liquid of linear molecules the semischematic mode-coupling model, previously intro- 
duced to describe the center of mass (COM) slow dynamics of a network-forming molecular liquid. 
We compare the theoretical predictions and numerical results from a molecular dynamics simulation, 
both for the time and the wave-vector dependence of the COM density-density correlation function. 
We discuss the relationship between the presented analysis and the results from an approximate so- 
lution of the equations from molecular mode-coupling theory [R. Schilling and T. Scheidsteger,P/i?/s. 
Rev. Em 2932 (1997)]. 
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I. INTRODUCTION 

The mode-coupling theory (MCT) |Q,E| has opened 
new perspectives in the theoretical understanding of 
the dynamic slowing down characteristic of supercooled 
glass-forming liquids MCT, originally developed to 
describe the structural relaxation in simple liquids, i.e. in 
liquids composed by particles interacting via spherically 
symmetric intermolecular potentials, also offers a coher- 
ent picture of the slow dynamics in molecular glass form- 
ing liquids, composed of asymmetric molecules. The abil- 
ity of MCT to model the onset of slow dynamics in molec- 
ular liquids has prompted the need to extend the theory 
to fully take into account the angular degrees of free- 
dom. The extension of MCT to a solute linear molecule 
in a solvent of spherical particles [Q and to molecular 
liquids [|||] , which we refer to in the following as molec- 
ular MCT (MMCT), has been recently achieved. The 
center of mass (COM) density-density correlation func- 
tion, which in MCT is the only relevant correlation func- 
tion, in MMCT becomes coupled to an infinite hierarchy 
of rotational correlation functions, arising from the ex- 
pansion of the angular degrees of freedom in spherical 
harmonics. The MMCT equations for COM and angular 
correlators have been solved until now for systems of lin- 
ear molecules under specific approximations Work 
is currently underway to improve the approximations for 
dumbbells and to calculate a solution for the general case 
of molecules of arbitrary shape. 

Recently, some of us proposed a parameterization of 
the role of the rotational degrees of freedom and their 
effective coupling with the COM density ||. This ap- 
proach, which provides a solvable set of coupled equa- 
tions for the slow dynamics of the COM density-density 
correlation function, has been named semischematic be- 
cause it retains all the ^-dependence of the COM cor- 
relators but condenses the coupling between COM and 



angular correlation functions into a single q-indepcndcnt 
parameter xr- A detailed comparison between the theo- 
retical predictions of the model and results from a molec- 
ular dynamics (MD) simulation has been performed for 
a network forming liquid j8|, finding an excellent agree- 
ment up to a cut-off q- vector value where the microscopic 
geometric details become dominant. 

In the long term, the exact solutions of the MMCT 
equations are to be preferred because they do not re- 
quire an ad hoc estimate of the translational-rotational 
coupling and predict also the behavior of the angular cor- 
relators. On the other hand, in the short term the semis- 
chematic equations are rather appealing because of their 
simplicity, for the minimal amount of input information 
required and because the complete time dependent solu- 
tion can be achieved with present day computational fa- 
cilities. Also, once the ability of the semischematic model 
to describe the time evolution of the COM correlation 
functions has been assessed by detailed comparison with 
MD simulations, comparison with experimental data is 
foreseeable, again due to the limited need for structural 
input. For this reason, in this Article we solve the semis- 
chematic model for a system of linear dumbbells inter- 
acting through a Lennard- Jones (LJ) potential and com- 
pare the predictions of the model with the correspond- 
ing quantities evaluated from long MD simulations of the 
same liquid ||. The choice of a liquid of dumbbells for 
which the MMCT equations have been previously solved 
approximately @ allows at the same time a comparative 
discussion of the two theoretical methods. 

II. THEORY 

The semischematic model is defined by introducing in 
the ideal MCT equations for simple liquids an effective 
coupling parameter \R which models the caging effect 
of the molecular rotational motion on the COM dynam- 
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ics, i.e. the slowing down of the COM relaxation intro- 
duced by the angular degrees of freedom. The resulting 
system of integro-differential equations [|| describes the 
time evolution of the normalized density-density COM 
correlation functions 4> q {t) — S q (t)/S q , S q (t) being the 
dynamic structure factor 
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and S q — S q (0) the static structure factor. The unknown 
Xr is fixed once and for all by requiring that the ideal 
glass transition temperature in the model T^ CT coin- 
cides with the temperature calculated form the analysis 
of experimental or MD data. The imposed equality of 
the theoretical and MD temperatures is very important 
because it allows to compare the numerical and the MD 
data at the same reference T, i.e. with the same structure 
factor and the same thermal broadening |ic|j . A larger 
value of xr corresponds to a stronger slowing down of 
the COM relaxation due to the interaction with the ro- 
tational motion. If xr — 1 the semischematic model 
coincides with the standard MCT 

The physics described by the semischematic model is 
the same as for the ideal MCT. It predicts that when 
T < T c MC " r , (f> q (t) does not relax to zero at long times 
and the COM dynamics is frozen. Above T C MCT the cor- 
relator 4> q (t) decays with a typical two step relaxation 
process, characterized by the fast decay to a plateau value 
(the non ergodicity parameter) followed by a slow relax- 
ation to zero which gets slower and slower as T C MCT is 
approached (a- relaxation) . 

The semischematic model describes the dynamic evolu- 
tion of (j) q (t) in the time region where the slow dynamics 
becomes dominant (a-region) by the system of coupled 
equations 
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4> q (t) = m q (t) -z / ds m q (t - s)4> q (s) 

dt J a 
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where the time variable i is defined in terms of a char- 



acteristic time scale which diverges at T t 
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memory function m q is a quadratic functional of the cor- 
relations 4> q {t) themselves 
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and its increase on cooling is responsible for the slowing 
down of the relaxation process. The parameter xr. enters 
in Eq. | as a ^-independent multiplicative factor, thus 
increasing the strength of the COM memory function. 
The vertices in Eq. |^ are defined as 



V{q,k)=S q S k S w _ % ^ 



q- k(l - S^ 1 ) + q- (q- fc)(l - S^) 
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and are functions of the COM S q and inversely propor- 
tional to the density. 

The infinite time limit f q of <j> q (t) can be calculated 
solving self-consistently on a discrete set of q values the 
coupled integral equations 
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using as input the number density n and S q which can 
be calculated from the MD data or measured experimen- 
tally. 

Having chosen xr properly, at the ideal glass transition 
temperature the solution f q of Eq. || jumps discontinu- 
ously from zero to a nonzero value, which defines the non 
ergodicity parameter /| 0J^] . In what follows we will ne- 
glect the upper index in the non ergodicity parameter. 

In the supercooled liquid phase, after the plateau f q , 
4> q (t) follows an initial power law decay (von Schweidler 
law), ruled by a q independent scaling exponent b, fol- 
lowed by a stretched exponential relaxation (Kohlrausch- 
William- Watts law): 



6,(i)~ A* exp 
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The range of validity of the von Schweidler law is strongly 
q dependent p"l) and, therefore, it is worthwhile to con- 
sider also the second order corrections: 
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The exponent b can be calculated solving 
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where T is the Euler gamma function and the exponent 
parameter X is defined by 
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In Eq. [9] e c and e c are the right and left eigenvector cor- 
responding to the maximum eigenvalue of the stability 
matrix 
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evaluated at the critical point. The critical amplitude 
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while /ig 2 ' can be calculated as explained in Ref. |I I j 

In Sec. IV we compare the predictions of the semis- 
chematic model for the slow relaxation in a liquid of linear 
molecules with the corresponding predictions of MMCT. 
To facilitate the reading of the present Article we present 
here a brief outline of the MMCT for linear molecules. 
We refer the interested reader to Refs. §-0]. 

According to MMCT the relevant variables for the de- 
scription of the slow dynamics of supercooled molecular 
glass forming liquids are the generalized correlation func- 
tions 

S${q,t) = 1 < Pi m (q,t)*pv m (q,0) > (12) 
where the density functions are defined by 



S]J l (q) up to I — 2 as evaluated from the set of MD data 
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In Eq. the sum runs over the N molecules of the liquid, 
Xj is the COM position of the jih molecule and Yi m (Clj) 
are the spherical harmonics for its Euler's angles. The 
index I ranges over the whole set of nonnegative integer 
numbers, while m € [—/,/]• For I = I' = the correla- 
tion function defined in Eq. [l^ coincides with the COM 
correlator S q studied by the usual MCT and the semis- 
chematic model. Due to the lack of rotational invariance 
of the molecules the correlation functions depend both 
on modulus and orientation of the wave vector. In Eq. 
|l2| we have chosen a reference frame where q points in 
the direction of the z axes, which offers the advantage 
of diagonality of the correlators with respect to m. The 
evolution equations for the correlators in Eq. 
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are a 

generalization of Eq. || for a single q dependent corre- 
lator. The slowing down of the relaxation is ruled by 
a infinite set of memory functions [MjJ}(q, t)] aa which 
are quadratic functionals of the whole set of correlators 
(|l2|). The extra indices a, a' £ {T, R} are related to pro- 
jection operations on the longitudinal translational (T) 
and rotational (R) currents. For the aims of this Article 
we underline that in MMCT the time evolution of each 
correlator is coupled, through the memory functions, to 
every other correlator. This means that, due to the de- 
pendence of [MQ (q, t)] TT on both the translational and 
rotational correlators, the dynamics of the COM corre- 
lator SQ (q,t) is affected by the time evolution of every 
angular correlation function S™, (q,t) with I ^ 0. 

With the cut-off I < 2, the non ergodicity parame- 
ters f 'u l {q) — lim^oo S" l l (q,t)/S'^(q) have been calcu- 
lated for the liquid of LJ dumbbells in the approxima- 
tion of diagonality in I, i.e. S™,(q,t) ps Su> S™(q,t) and 
[MJJ}(q,t)] aa ' re S w [MJJ l (q,t)] aa ' §. We stress that this 
diagonality was also demanded for the static correlators 
SJu(q), in order to keep the MMCT-equations as simple 
as possible. The input of the calculations are the number 
density of molecules and all the static structure factors 



In Fig. |] we show the results (from Ref. |^|) for 
the COM non ergodicity parameter /q at two different 
temperatures (short and long dashed lines). The higher 
one, T t MMCT = 0.383, corresponds to the MMCT theo- 
retical freezing temperature for the translational degrees 
of freedom, i.e. at T t MMCT the COM / ° jumps to a non 
zero value while the rotational dynamics is still in a liq- 
uid phase. At the lower temperature T r MMCT = 0.310 all 
the other non ergodicity parameters f™ (I ^ 0) become 
nonzero, i.e. T r MMCT is the MMCT theoretical transition 
temperature for the rotational dynamics. This splitting 
of critical temperatures, not observed in the simulation, 
is an artifact of the diagonalization approximation and 
is not observed when the approximation I ^ I' is waived 



III. MODEL AND SIMULATION DATA 

The model under investigation is an one-component 
system containing N — 500 rigid diatomic molecules. 
Each molecule consists of two atoms, labeled A and B, 
separated by a distance d. The interaction between two 
molecules is build up by pair interactions between the 
atoms, which are due to the LJ potentials 



V a p(r) = 4e a/3 



a,/3 e {A,B} 
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with LJ parameters eaa — £ab — 1-0, ess = 0.8 and 
o aa = cab = 1-0, &BB = 0.95, i.e. 6,4^4 was chosen 
as the unit of energy and temperature (fcs = 1) and 

o aa as the unit of length. The unit of time is then 

i 

((a AA m)/(4:8eAA)) 2 where m is the mass of an atom 
which is chosen to be equal for both types of atoms. The 
slight head-tail asymmetry of the dumbbell assures, to- 
gether with the choice of d = 0.5 as inter-atomic dis- 
tance, a good coupling between translational and rota- 
tional motion on the one hand and avoids crystallization 
into a liquid crystalline phase and the intersection of two 
dumbbells on the other. 

After equilibrating the system in the (N,p,T)-ensemble 
for times which exceeded the relaxation times of the sys- 
tem even at the lowest temperature, the production runs 
were carried out in the microcanonical ensemble. To im- 
prove the statistics the data for each temperature were 
averaged over at least 8 independent runs. Further de- 
tails about the simulation can be found in the references 
[^| from which we take part of the data to be compared 
with the theoretical results. 

The critical temperature for the COM which was de- 
termined from the simulation by fitting the a-relaxation- 
time and the diffusion constant D with power laws 
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roc(T-T c )-^ , Dot{T-T c y (15) 

according to MCT predictions, was found to be T^' ID = 
0.477. The numerical values for the non ergodicity pa- 
rameter f q , the critical amplitude hq 1 ^ and the second 

order correction h q 2 ^ at the critical temperature T MD 
0.477 are also taken from reference 0. They were evalu- 
ated by fitting the decay from the plateau in <fi q (t) with 
the von Schweidler law Eq. [?], including the second order 
correction. From the same procedure the critical expo- 
nent b = 0.55 is also obtained and, via Eq.||, one gets 
for the exponent parameter A the result A = 0.76. Fur- 
thermore we have examined the time dependence of the 
MD correlator <p q (t) in the a-region by evaluating the 
amplitudes Aff , stretching exponents 0^ and relaxation 
times Tq of a Kohlrausch- Williams- Watts fit (Eq. |6|). 

IV. RESULTS 

We solve for the liquid of LJ dumbbells the semis- 
chematic equations introduced in Sec. [H] using as in- 
puts of the calculation the static structure factor S q as 
obtained from the simulation (Figs. ||||) and the COM 
number density n = 0.719 at T MD = 0.477. By solv- 
ing Eq.|^ on a grid of 300 equispaced q vectors, we find 
that the condition T C MCT = T MD fixes the value of xr to 
1.17, which suggests that the coupling between COM and 
angular degrees of freedom in LJ dumbbells increases the 
COM memory function about 20%. If compared to the 
value xr = 1-93 found for SPC/E water |Q, this result 
highlights the weaker hindering effect of the rotational 
motion in a liquid of LJ linear molecules with respect to 
the strong one observed in an hydrogen-bonded network- 
forming liquid. In water the highly energetic hydrogen 
bonds build up the network structure which is responsi- 
ble for the caging of the molecules in the glass phase. The 
motion of the COM of the molecules is completely slaved 
to the breaking and reforming of the hydrogen bonds, i.e. 
it is definitely correlated to the angular dynamics. 

The solution of Eq|| for xr = 1-17 as a function of 
the wave-vector q is shown in Fig. [I] together with the 
COM non ergodicity parameters as calculated from the 
MD data. The theoretical f q oscillates in phase with the 
MD data, but underestimates the amplitude, especially 
at large q vectors. The shoulder around q = 3, which may 
be attributed to the rotational-translational coupling, as 
the orientational correlator Si^q) has a maximum at 
this q, is also underestimated. This notwithstanding, the 
semischematic f q captures the q-dependence of the non 
ergodicity parameter as calculated from the MD data. 

The simplicity of the semischematic equations allows 
us to study, besides f q , also the complete time relax- 
ation of <j) q (t). According to the theoretical predictions 
outlined in Sec. O we calculate the critical amplitudes 



h q and h q ' and the exponent b which rule the early 
a-relaxation behavior (Eq. [?]) . 

In Fig. ||, we present a comparison between the the- 
oretical predictions of the semischematic model for the 
critical amplitudes hq and h q 2 ^ and the same quantities 
as calculated by fitting Eq.f?] to the density-density cor- 
relation functions, i.e. with a quadratic fit in t b . The 
fitting coefficients h q /r b and h q 2 ^ /r 2b can be compared 
with the theoretical critical amplitudes after fixing, once 
and for all, the q-independent time scale r, introduced 
by the scale invariance of Eq. ||. The agreement is satis- 
factory in a wide range of q values. The theoretical value 
for the exponent parameter, as calculated from Eq. 
is A = 0.63 to be compared with A = 0.76 as obtained 
by the simulation, i.e. the difference between the two is 
about 15%. This yields a theoretical critical exponent 
b = 0.75, while the exponent calculated from MD data is 
b = 0.55. In the case of SPC/E water the theoretical and 
MD values of b coincide within the numeric error. Such 
finding is consistent with the remarkable agreement of the 
non ergodicity parameters over both the relevant peaks 
of the structure factor. 

The comparison between MD data and prediction of 
the semischematic model can be extended to the long 
time region. We solve the complete dynamic set of equa- 
tions ^ in the whole q range We fit the stretched 
exponential law (Eq. ^) to the long time relaxation (late 
a-region) and we compare the obtained amplitude A^ , 
relaxation time and stretching exponent f3jf with the 
corresponding quantities as calculated by fitting the MD 
relaxation. The comparison of the complete q depen- 
dence of Aq , Tq and (3^ is shown in Fig. ||. The 
theoretical and numerical relaxation times are in per- 
fect agreement. Less satisfactory, as in the case of water, 
is the theoretical prediction for f3q for which the the- 
ory provides the correct qualitative q-dependence, but 
failing in amplitude up to 30%. The error in the values 
of the stretching exponents is expected on the basis of 
the drastic simplification adopted in the semischematic 
approach, which condensate the coupling between the in- 
finite set of angular correlator and the COM correlator. 
Indeed, in phenomena in which the decay of correlation 
results from the sum of several independent relaxation 
processes, a smaller j3 indicates a wider distribution of 
relaxation times (lfj. As expected A^ has the same be- 
havior as f q both for theory and simulation. 

The choice of the dumbbell liquid is particularly in- 
teresting because it allows a comparison between the 
theoretical predictions of the semischematic model and 
those of MMCT which provides a deeper understanding 
of the basic approximation in the model, i.e. the assump- 
tion that the coupling between the COM and angular 
degrees of freedom can be quantified in a multiplicative 
q-independent factor xr- The comparison requires a cer- 
tain degree of care, because of the difference in the theo- 
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retical estimate of T MMCT , rather different from T^ D . 
Indeed, while in the case of the semischematic model the 
MD and the theoretical f q are evaluated at the same 
temperature, in the case of MMCT, the non ergodicity 
parameters are calculated using as inputs the structure 
factors from the simulations but evaluated at tempera- 
tures different from T C M£ \ Fig. reports the COM 
MMCT non ergodicity parameters obtained in the ap- 
proximations of diagonality and I < 2 for the freezing 
of the COM dynamics only, and for the freezing of both 
COM and angular degrees of freedom. 

A better insight into the comparison between the dif- 
ferent theoretical approaches can be performed studying 
the non ergodicity parameter as predicted by the semis- 
chematic model at T^' IMCT . Keeping fixed once and for 
all the coupling xr — 1-17, we can solve Eq. || with the 
memory function Eq.[| With this choice of temperature 
the COM dynamics is in a deep glassy phase, both in the 
semischematic and MMCT descriptions. Thus, being in 
a non ergodic phase, the COM dynamic structure fac- 
tor has a finite long-time limit, which is shown in Fig. 
|l|. The semischematic (long dashed line) and MMCT 
(dotted line) predictions are in perfect agreement. This 
result can be enlightened by the comparison of the long 
time limits of the semischematic memory function m q 
and the COM MMCT memory function [M ° (<?)] TT at 
T MMCT ( Fig ^ We recall that the angu i ar correlators 

contribute to [AlQ Q (q)] TT as well as the COM correlator, 
while the semischematic m q is ruled only by the COM 
S q and the value of XR- Thus, the semischematic model 
gives rise to the complete functional dependence of the 
COM memory function on the rotational relaxation in a 
very simple way, i.e. taking into account only the func- 
tional dependence on the COM correlator and summariz- 
ing all the remaining coupling in the effective xr- This 
means that the COM dynamics predicted by the semis- 
chematic model almost coincides with the corresponding 
predictions of MMCT in the diagonalization approxima- 
tion. Furthermore the model has the advantage that the 
complete time evolution of the COM correlator can be 
calculated with current computational resources. 



that the q-dependence of S q can be reduced to a sin- 
gle representative qo vector, i.e. S q oc S(q — qo) fL8[ . 
In this approximation, which is named schematic model, 
the g-dependence is abandoned in favor of an exact de- 
scription of one or two representative correlators. In the 
same spirit, the semischematic model provides a method 
for studying the complete g-dependence of the COM dy- 
namics, neglecting the angular degrees of freedom which 
are condensed in a single g-independent parameter xr- 

Notwithstanding the drastic approximation intrinsic in 
the proposed approach, we have shown in this Article 
that the model captures the essential ingredients of the 
q dependent static and dynamic features of the COM a- 
relaxation. The predicted non ergodicity parameter, re- 
laxation time and stretching exponent oscillate in phase, 
the critical amplitudes out of phase, with the COM struc- 
ture factor, i.e. the model predicts the same qualitative 
behavior observed in the simulation. The quantitative 
agreement with the MD data for f q , h q X \ h q 2 ^ and T q is 
satisfactory, especially for q vectors close to the maxi- 
mum of S q , while the stretching exponent [3k is overesti- 
mated. This reflects the major weakness of the approach 
but at the same time clearly indicates the role played by 
the angular degrees of freedom in controlling the dynam- 
ical evolution of the center of mass. If the present ob- 
servations are discussed together with the semischematic 
analysis of the dynamics of SPC/E water, a model which 
mimics a liquid of strong directional hydrogen bonds, it 
gets obvious that the value of xr is a measure for the 
strength of the roto-translational coupling. 

The comparison between the semischematic model and 
the MMCT approximate solutions for the non ergodicity 
parameter and the memory function also supports the 
validity of the assumption of a g-independent xr- More- 
over, the coincidence of the theoretical critical T and the 
MD one, opens the way to precise comparisons between 
theory and experiments. Differently from the MMCT, 
the only required input in the model is the COM S q , a 
quantity which can be experimentally measured by suit- 
able designed neutron or x-rays scattering experiments 



V. CONCLUSIONS 

A theoretical (but relatively simple) description of the 
COM dynamics in a molecular supercooled liquid is a 
relevant task which is recently focusing a lot of efforts. 
The situation recalls the early day of the MCT for sim- 
ple liquids, when the dynamical equations were known 
but the exact numerical solutions were too difficult to 
handle. In that situation, simple approximations giv- 
ing rise to a solvable set of equations were proposed and 
carefully studied. This class of approximations, which is 
still extensively used to interpret in a simple way exper- 
imental results 16 17j], arises from the basic assumption 
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FIG. 1. COM non ergodicity parameter f q as calculated by solving Eq. q (solid line) and as evaluated from the MD 
simulations j^] (symbols). For comparison also the COM non ergodicity parameters as predicted by MMCT are reported at 
two different temperatures (dashed lines). The dotted line shows f q as predicted by the semischematic model at T — 0.310. 
The unit of q is 2tt(t a ^ a while the unit of temperature is £aa (setting Kb = 1). The MMCT curves are from Ref. ||. 
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FIG. 2. Critical amplitudes and in the von Schweidler law Eq.^ as predicted by the semischematic model (solid 
lines) and as calculated from the MD data (symbols). The static structure factor S q used as input of the theoretical calculation 
is also shown. The unit of q is 2iva^. 



8 




FIG. 3. Amplitude A^, stretching exponent f3^ and relaxation time in the Kohlrausch- William- Watts law Eq. g. Lines 
are obtained fitting Eq. pi to the exact time dependent solution of Eq. ti while symbols are evaluated fitting the same law to 
the MD data. The MD time unit for is f0 5 ((cr^ A m)/(48eAA)) 1 ^ 2 while the MCT relaxation times are arbitrarily scaled. 
The unit of q is 2-na^\. The static structure factor is shown as reference. 
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FIG. 4. COM memory function as predicted by the semischematic model Eq. |3| (solid line) and by MMCT (dashed line) at 
rpMMCT _ o 310 'pjjg un [i of q is 2na^ A while the unit of temperature is eaa (setting Kb = !)• 
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